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ABSTRACT 

The Bondi-Hoyle formula gives the approximate accretion rate onto a point particle accreting from a uniform 
medium. However, in many situations accretion onto point particles occurs from media that are turbulent rather 
than uniform. In this paper, we give an approximate solution to the problem of a point particle accreting from 
an ambient medium of supersonic ally turbulent gas. Accretion in such media is bimodal, at some points re- 
sembling classical Bondi-Hoyle flow, and in other cases being closer to the vorticity-dominated accretion flows 
recently studied by Krumholz, McKee, & Klein. Based on this observation, we develop a theoretical prediction 
for the accretion rate, and confirm that our predictions are highly consistent with the results of numerical simu- 
lations. The distribution of accretion rates is lognormal, and the mean accretion rate in supersonically turbulent 
gas can be substantially enhanced above the value that would be predicted by a naive application of the Bondi- 
Hoyle formula. However, it can also be suppressed by the vorticity, just as Krumholz, McKee, & Klein found 
for non-supersonic vorticity-dominated flows. Magnetic fields, which we have not included in these models, 
may further inhibit accretion. Our results have significant implications for a number astrophysical problems, 
ranging from star formation to the black holes in galactic centers. In particular, there are likely to be significant 
errors in results that assume that accretion from turbulent media occurs at the unmodified Bondi-Hoyle rate, or 
that are based on simulations that do not resolve the Bondi-Hoyle radius of accreting objects. 
Subject headings: Accretion, accretion disks — black hole physics — hydrodynamics — stars: formation — 
turbulence 



1. INTRODUCTION 

Accretion of gas from a background medium onto a small 
particle occurs throughout astrophysics. Examples range 
from protostars accreting from their natal clouds to black 
holes accreting interstellar gas during galaxy mergers. In 
such cases, one wishes to know the rate of mass accre- 
tion onto the point-like object. Bondi, H oyle, and Lyttle- 
ton (iHovle & LvttletonlfT939t lBondilll952l) derived the clas- 
sic solution to this problem when the background gas is uni- 
form and either stationary or moving with constant veloc- 
ity relative to the accreting object. More recent numerical 
work has demonstrated that their treatment, while missing 
some instabilities that cause the accretion rate to be time- 
dependent when the background gas is moving rapidly rela- 
tive to the accretor, gives essentially the correct accretion rate 
(iFrvxell & Taamlll988tlRuffertlll994lRuffert & Arnettlll994t 
iRufferl 19971119991: lFoglizzo"& Ruffert 19991) The a ccretion 
rate is reasonably well-fit by (Ruffert & Arnett 1994) 
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where p^, c^, and M. are the density, sound speed, and 
Mach number of the flow far from the accreting object, and 
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A is a factor of order unity that depends on M. and the equa- 
tion of state of the gas. For M. = (Bondi accretion) in an 
isothermal medium, A = exp(1.5)/4 w 1.1, and we adopt this 
value in our analysis below. 

However, many sources of gas accretion are supersonically 
turbulent, and thus are not easily characterized by a constant 
background density or velocity. A seed protostar accreting 
gas to which it is not bound in a process of "c ompetitive 
accretion" (e.g. Bonnell et al. 1997; Klessen & Burkert 2000, 
1200 11 iBonnell et al JI200 1 albl 120041) inside a tur bulent molec- 
ular co re is an example of such a phenomenon. Pado an et alJ 
(2005) suggest that the observed accretion rates onto T Tauri 
stars may be the result of accretion from the remnant molecu- 
lar gas of the star-forming region. 

Despite the ubiquity of the phenomenon, theoretical and nu- 
merical work to date has not clearly shown how to extend the 
Bondi-Hoyle-Lyttleton result to the case of a turbulent back- 
ground medium. As a zeroth order estimate, one might simply 
use the mean density and gas v elocity dispersion of the turbu- 
lent medium in the Bondi ( 119521) formula in place of p^ and 
Coo. This would produce a zeroth-order estimate of the accre- 
tion rate 
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where p is the mean density, A4q is the Mach number of the 
turbulent flow, c s is the sound speed, and M is the mass of the 
accreting object. Of course this approximation is only valid 
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for M £ 1 ■ fradoan et alJ J2005h suggest a somewhat more so- 
phisticated approach in which one computes the accretion rate 
by applying the Bondi-Hoyle formula at every point within a 
turbulent medium and then taking the volume average. 

However, numerical and analytic studies have shown that 
even a small amount of vorticity in the a ccreting medium can 



substantially c hange the accretion rat e (Scarke & Sh ulll98C 
ISoarkd 1 19821 lAbramowicz & Zurekl 119811 iFrvxeli & Taaml 
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1999; Igumenshchev & Abramowicz 
1 20001 llgum enshche v et alJ 120001 iProga & BegelmaiJ 12003 
Krumholz et al. 2005| in regions where the overall flow ve- 
locity is small, and in a supersonically turbulent medium the 
bar oclinic i nstability is like ly to generate significant vortic- 
ity. iRrumholz et alJ ( 120051) use simulations and analytic ar- 
guments to derive a formula analogous to the Bondi-Hoyle 
formula accretion in a medium with vorticity, and application 
of this formula to turbulent media suggests that the vorticity 
of a turbulent gas may be at least as important as its veloc- 
ity in inhibiting accretion. Magnetic fields, which we do not 
include here, may further reduce the accretion rate. 

In this paper we derive an estimate for the accretion rate 
onto a point particle in a turbulent medium, thereby extend- 
ing the Bondi-Hoyle-Lyttleton solution to this case. We limit 
ourselves to the case of an accretor substantially smaller than 
the medium from which it is accreting, where the self-gravity 
of the medium is negligible in comparison to the gravity 
of the accretor, and where the accreting medium is dom- 
inated by supersonic, isothermal turbulence. Stars accret- 
ing in molecular clouds satisfy these conditions due to the 
short times the gas requires to reach radiative equilibrium 
( Vazquez-Sema deni et al.ll2000l) . and below we discuss other 
astrophysical situations to which our resuls apply. In §|2]we 
propose a simple method for determining the accretion rate in 
a turbulent medium as a function of the properties of the gas. 
In §|3] we describe numerical simulations we have conducted 
to test our theoretical model, and show that it provides an ex- 
cellent fit. In contrast, earlier proposed models give far less 
accurate predictions. In §|4]we use our model to predict the 
rate of accretion onto point particles in a turbulent medium 
as a function of a few simple properties of the accreting gas, 
and in §[5] we discuss the implications and limitations of our 
approach. Finally, in §[6] we summarize and present our con- 
clusions. 

2. A SIMPLE MODEL FOR ACCRETION IN A TURBULENT FLOW 

In this section, we describe a simple ansatz to relate the 
rate of accretion onto a point particle in a turbulent gas to the 
properties of the gas. Since a turbulent medium has fluctua- 
tions in gas properties in space and time, it is most convenient 
to characterize the accretion rate for a particle placed within 
it using a cumulative probability distribution function (PDF), 
which specifies the probability P^(< M) that a particle placed 
at a random position within the turbulent gas will accrete at 
a rate less than M. We wish to predict this function in terms 
of the properties of the turbulent gas. (Throughout this pa- 
per, we use P to refer to cumulative distribution functions and 
dp to refer to the corresponding differential distribution, or 
probability density, functions.) 

Our ansatz is that we can roughly divide accretion in a tur- 
bulent medium into two modes. In regions with relatively 
high velocities and small vorticities, the vorticity should have 
little effect on the flow pattern, and flow shou ld resemble or- 
dinary Bondi-Hoyle accretion. For example, Ruffert (1997, 
1999) show that in simulations with relatively little vorticity, 



the overall accretion rate and flow pattern are quite similar to 
what one finds with no vorticity. In regions where the flow 
velocity is relatively small and the vorticity is relatively large, 
the small overall velocity of the gas relative to the particle 
should have little effect. Instead, the vorticity should domi- 
nate the flow, and the flo w pattern and accretio n rate should 
resemble those found by Krumholz et al] (120051) . 

Which mode occurs for a given particle will be determined 
by which one would produce a lower accretion rate. Thus, our 
suggested procedure for estimating P^{< M) in terms of the 
properties of the turbulent gas is to compute the function 

M turb (x) « [M BH (x)- 2 +M tJ (x)- 2 ]" 1/2 , (3) 

at every point x within the flow. Here, Mbh is to be com- 
puted with the Bondi-Hoyle formula (equation 0, using the 
density p(x) and velocity v(x) at x for p^ and Vqq, and us- 
ing the constant isothermal sound speed c s . The quantity M w 
is the accretion rate in a vorticity-dominated medium. It is 
a function of the den sity and vorticity u> = |V x v|, given by 
JKrumholz et all2005l) 
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is the dimensionless vorticity, 
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is the Bondi radius of the accreting object , and the func- 
tion f (u)*) is given in terms of an integral in Krumho lz et alJ 
(2005). For our purposes it is convenient to approximate the 
integral by 

/<«.)« J^, (7) 

which is accurate to better than 12% for lq* < 10 4 . Equation 
Q simply computes the harmonic mean of the squared accre- 
tion rates, effectively minimizing between them. 

We could test our ansatz by simulating regions with a range 
of vorticities and velocities, but that parameter space is rel- 
atively large and exploring it would be time-consuming. In- 
stead, we can test our model directly against a turbulent re- 
gion. If such a region has volume V, then we predict that 

Pm« M) = ^ / H[M-M tmh (x)] dV, (8) 

where H(x) is the Heaviside step function, which is unity for 
x > and zero for x < 0. In contrast, the model proposed by 
IPadoan et alJ (120051) suggest using Mbh (x) instead of M tur b (x) . 
Note that this equation implicitly assumes that the we can es- 
timate the accretion rate onto a particle by looking at the dis- 
tribution of gas properties in space rather than the distribu- 
tion of gas properties in time at the particle's location, which 
might seem more intuitive. However, since turbulent media 
are roughly homogenous in time and space (at least over time 
scales shorter than the time it takes the turbulence to decay 
and on length scales smaller than the outer scale of the tur- 
bulence) the space and time distributions should be approxi- 
mately the same. We can test this homogeneity approximation 
by comparing our predictions to simulations, as we do in §[3] 
Because we are assuming that the spatial and temporal dis- 
tribution of accretion rates are the same, for much of what 
follows we will not explicitly distinguish between the two. 
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3. SIMULATIONS 

3.1. Simulation Methodology 

To test our model, we simulate accretion onto point parti- 
cles in a turbulent medium. Our calculations use our three- 
dimensional adaptive mesh refinment (AMR) code to solve 
the Euler equations of compressible gas dynamics 
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where p is the density, v is the vector velocity, P is the thermal 
pressure (equal to pc 2 s since we adopt an isothermal equation 
of state), e is the total non-gravitational energy per unit mass, 
and 4> is the gravitational potential. The code solves these 
equations using a conservative high-order Godunov scheme 
with an optimized approximate Riemann solver (Toro 1997), 
and its implement ation i s described in detail in Truelov e et alJ 
( 1998) and lKleir] ill 9991) . The algorithm is second-order accu- 
rate in both space and time for smooth flows, and it provides 
robust treatment of shocks and discontinuities. We adopt a 
near-isothermal equation of state, so that the ratio of specific 
heats of the gas is 7 = 1.001. 

Although the code is capable of solving the Poisson equa- 
tion for the gravitational field <fi of the gas based on the density 
distribution, as discussed above we neglect the self-gravity of 
the gas and consider only the gravitational potential of the ac- 
creting particles. Thus, the potential is given by 
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where M,- and x, are the mass and position of particle 2. 
The particles themselves are Lagrangia n sink particles imple- 
mented using the algorithm of Krumho lz et alJ (|2004). They 
are capable of moving through the gas, accreting from it, and 
interacting gravitationally with the gas. Although ordinar- 
ily sink particle can also gravitationally interact with one an- 
other, in this simulation we neglect inter-particle gravitational 
forces. 

The entire code operates within the AMR framework 
(Berger & OligeJ 11984 iBerger & Collelal 119891: iBell et alJ 
1984). In AMR, one discretizes the problem domain onto a 
base, coarse grid, which we call level 0. The code then dy- 
namically creates finer levels L = 1,2, . . .n within this based 
on user-specified criteria. The cells on level L have cell spac- 
ings that are a factor of / smaller than those on the next 
coarser level, i.e. Axi = Ax^-i/f, where / is also user- 
specified. For the runs in this paper, we use / = 2. To take 
a time step, one advances level through a time Afo, and then 
advances level 1 through / steps of size At\ = Ato/f. How- 
ever, for each advance on level 1, one must advance level 2 
through / time steps of size A?2 = At\/f, and so forth to the 
finest level present. The flux across a boundary between level 
L and level L + 1 computed during one time step on level L 
may not match that computed over / time steps on level L+ 1 . 
For this reason, at the end of each set of / fine advances, the 
code performs a synchronization procedure at the boundaries 
between coarse and fine grids to ensure conservation of mass, 
momentum, and energy. 



3.2. Simulation Setup 

The initial conditions for our simulation consists of a box 
of gas with a uniform density p = 1 (we use non-dimensional 
units throughout) and sound speed c s = 1. The box has pe- 
riodic boundary conditions and extends from -1 to 1 in the 
x, y, and z directions. We impose an initial turbulent ve- 
locity field by genera ting a Gaussian-random grid of vectors 
( Dubinski et al. 1995). The grid has all its power at wavenum- 
bers k = 1-2, where for a wave of wavelength A, k = 4/ A. 
Thus, k = 1 corresponds to the largest wave that will fit in the 
box. We normalize the initial velocity so that the 3-D Mach 
number is .Mo = 40. Once the we have established the ini- 
tial conditions, we allow the turbulence to decay freely un- 
til the Mach number reaches Ado = 5.0. The decay allows 
the turbulence to cascade down to small scales and reach the 
commonly observed k~ 2 spectrum for supersonic turbulence, 
which is the natural consequence of Burger's equation. Dur- 
ing this phase of the simulation we use a resolution of N = 5 12 
cells per linear dimension with a fixed (non-adaptive) grid. 
We choose 40 as the initial Mach number because that allows 
us to let the turbulence decay through more than one e-folding 
before we measure its properties, ensuring that the flow has 
reached statistical equilibrium. 

Once the turbulence has decayed to A4q = 5.0, we insert 
sink particles into the simulation. The particles are arrayed 
in a uniformly spaced 4x4x4 grid, and have masses M = 
13/32 (in units where G = 1), giving them Bondi-Hoyle radii 
at A^ = 5.0 of 

M 1 .... 

rBH = - T72 = 77- ( 13 ) 

1+Ml 64 

Since the panicles are separated from their nearest neighbors 
by a distance of 1/2, they are separated by 32reH and should 
accrete completely independently of one another. When we 
introduce the particles, we refine the regions around them to 
guarantee that at every point in space r j Ax > 16, where r is 
the distance to the nearest particle and Ax is the grid spac- 
ing of the finest grid covering that point. We continue this 
refinement to a maximum resolution of N = 8192, giving a 
grid spacing of Ax = 1/4096 = tbh/64 on the finest level. We 
d iscus s issues of resolution and convergence in more detail in 

§EH 

3.3. Simulation Results 

After we have inserted the sink particles, we allow a sim- 
ulation to evolve until both the median and mean accretion 
rates (mean and median over particles, not over time) onto the 
particles becomes roughly constant, as shown in Figure^ We 
estimate by eye that the accretion rate has stabilized at time 
t > 4?bh, where 

fBH^xr- =0.00313. (14) 

The accretion rate in the plot is normlized to Mo = 0.0166, in 
our dimensionless units. 

As our simple model predicts, after the accretion rate has 
reached equilibrium the flow pattern around some of the par- 
ticles is similar to that of Bondi-Hoyle accretion with no vor- 
ticity, and the flow pattern around others is closer to the case 
of high vorticity. Figure [2] shows an example of the former. 
The velocity of the gas on the right side of the plot is rel- 
atively large and relatively uniform, indicating there is little 
vorticity. As the gas passes the sink particle, the sink parti- 
cle's gravity causes streamlines to converge and a shock to 



4 




0.04 



FIG. 1 . — Mean (upper line) and median (lower line) accretion rate versus 
time. The dotted vertical line shows the time when the accretion rate has 
reached equilibrium. 




FIG. 2. — Density (color) and velocity (arrows) around an accreting par- 
ticle. The origin is centered on the accreting particle, and the axes show 
position in code units and in units of reir- 

form. This leads to a clear Mach cone, as is normally seen in 
Bondi-Hoyle accretion (e.g. Krumholz et al. 2004, Figure 3). 
In contrast, Figure [3] shows a vorticity-dominated accretion 
flow. There is no Mach cone, but there is a dense, rotationally- 
supported torus of gas around the accreting particle. This is 
similar to the flow patterns that occur for accretion with no 
net relativ e velocity between the gas and the particle, just vor- 
ticity (e.g. Krumholz et al. 2005, Figure 5). 

We next compare the accretion rate predicted by our model 
to what we measure in the simulations. To measure the ac- 
cretion rate in the simulations, we wait until t > 4?bh, and 
compute the time-averaged accretion rate onto each particle 
thereafter. We then construct the cumulative distribution func- 
tion of the accretion rates onto the particles, which we wish to 
compare with P^(< M) as predicted by our models. For the 
model, we examine the simulation cube at the time when we 
insert the sink particles. We then compute Mbh (equation [2 
with A = exp( 1 .5)/4), M w (equation|4} and M tmb (equation[3} 
in every cell of the simulation. To measure vorticity, which 
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FIG. 3. — Density (color) and velocity (arrows) around an accreting par- 
ticle. The origin is centered on the accreting particle, and the axes show 
position in code units and in units of reH- 
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FIG. 4. — Cumulative (upper panel) and differential (lower panel) accre- 
tion rate PDFs vs. M/Mq. The curves shown are the measured value from 
the simulation (histogram) and the theoretical predictions using M t urb (solid 
line), Mbh (dashed line), and (dotted line). 

is not a primitive of our simulation, we take two-sided differ- 
ences over one cell at every point. From the grid of accretion 
rates, we compute P^{< M), using Mbh> M w , and M t m-b m 
equation (|SJl. We compare the predictions to the measured 
accretion rates in Figure 
Note that, since Jbh <C 1, the state of the simulation cube 
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on the large scale does not change much over the time after 
we insert the sink particles. For example, if we continue to let 
the box evolve without adding sink particles, after 4?bh, the 
Mach number only decreases from 5.0 to 4.7. Thus, there is 
no problem in applying our simple model at the time just be- 
fore we insert the sink particles. Also note that, at the time we 
measure vorticity, we have not yet inserted the sink particles 
and so there are no higher level AMR grids. When these grids 
first appear, the velocity is linearly interpolated on a cell-by- 
cell basis, so the vorticity in a refined cell is equal to that of its 
parent. However, as the simulation evolves, baroclinic insta- 
bility generates vorticity structure on more refined levels just 
as it does on the base level. By the time our accretion rates 
stabilize, the vorticity structure on the more refined levels has 
had several crossing-times to reach equilibrium. 

As Figure |4] shows, our ansatz function M tur b gives excel- 
lent agreement with the simulation results. In contrast, using 
either Mbh or M w in place of M tur b gives an extremely poor 
fit. We quantify this result by comparing the theoretical mod- 
els to the si mulations using a Kolmogorov-Smirnov (KS) test 
JPressetal]ll992l). The K S statistic for our model is 0.81. 
For the Padoan et al. ( 2005) proposal of using Mbh where we 
use M tul -b, it is 0.012, and using M w instead of M tur b gives 
2.7 x 10~ 10 . This indicates that our model is highly consis- 
tent with the simulation data, while the two other candidates 
we have checked are essentially ruled out. The difference is 
large enough to have substantial astrophysical significance: 
the accretion rate averaged over particles that we measure in 
our simulations is (M)/Mq = 0.73, and our theory predicts a 
mean of (M) /Mo = 0.90. (The difference is larger than one 
might expect because our 64 particles do not s ample the tail 
of the distribution very well.) In contrast, the Padoan et al. 
(2005) theory predicts (M) /Mo = 4.3, an ovestimate of nearly 
an order of magnitude. 

Our theoretical model is well-described by a fit to a lognor- 
mal PDF. We find a best fit of 

d PM ^n^^J [ln(M/M„)-l.l] 2 l 
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0.26 exp < 
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Figure |5] shows our fit in comparison to the simulation result. 

A final question to consider is whether we might be un- 
derestimating the accretion rate because we only count mass 
that falls onto the particle, not mass that accumulates in the 
rotationally-supported torus or other bound structures around 
the particle. One might worry that these structures are acting 
as mass reservoirs, and that the accretion rate might rise at 
some point in the future as mass accretes fro m them on a vis- 
cous ti me scale. A priori this seems unlikely: iKrumholz et alJ 
J2005I) simulate accretion with vorticity for 200 Bondi times 
and find that the rotationally-supported toruses do not act as 
reservoirs that raise the accretion rate at later times. Instead, 
the pattern is that the accretion rate starts high and gradually 
declines to its equilibrium value, while the toruses reach a 
steady mass, so that the ratio of torus mass to accreted mass 
declines steadily with time. Real accreting objects are likely 
to be in this limit where the torus mass is negligible compared 
to the total accreted mass. For example, the Bondi-Hoyle time 
for a protostar accr eting in from a turbulent molecular clump 
is, as we show in § 15.31 only ~ 10 yr, vastly shorter than the 
time for which the protostar accretes. 

We unfortunately cannot run our simulations for hundreds 
of Bondi-Hoyle times due to limitations of computational re- 
sources, and because over such long time scales the particles 
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FIG. 5. — Numerical (solid line) and best-fit theoretical model (dashed 
line) accretion rate probability distributions vs. M/Mq- 



would likely inte rfere with one an other. Thus, we cannot do 
a direct test as in IKrumholz et alJ J2005I) . However, to elim- 
inate the possibility that we are missing substantial accretion 
by neglecting the gas in bound structures, we compute all the 
mass within a distance tbh of each sink particle that has a po- 
tential energy larger than its kinetic energy, i.e. that is bound 
to that sink particle. We do this at 4^bh? when the particle 
accretion rate has reached equilibrium, and at the final time 
in the simulation, roughly 4.4?bh- We then compute the dif- 
ference between the bound mass at 4?bh and at the final time. 
The median change in the mass of bound structures, divided 
by the median change in mass accreted onto particles over the 
same period, is 0.27. Thus, even if we were to assume that 
bound structures were reservoirs of mass that grow with time, 
the extra accretion onto them would only represent a ~ 30% 
correction to the accretion rate we have computed. If we add 
the accretion rate onto bound structures to the accretion rate 
onto particles, and compare that to our theoretical model us- 
ing a KS test, we find a statistic of 0. 15. While this is not quite 
as good as our value obtained using just the particles, this is 
still a reasonably good fit. 

Note, however, that including the mass in bound struc- 
tures only makes sense if they are truly accreting. Otherwise, 
adding their mass to the mass accreted by the particles is sim- 
ply adding noise. It seems likely that this is the case, because 
the bound mass fluctuates strongly, and can be quite small 
compared to the accreted mass: for 16 of our 64 particles, the 
bound mass is less than 10% of the accreted mass at some 
point during the simulation, for 32 of the particles it is less 
than 25% of the accreted mass, and for 58 of them it is less 
than the accreted mass. Since the amount of mass that re- 
mains in the toruses over long times is small compared to the 
accreted mass, the toruses cannot be acting as mass reservoirs. 
Instead, the mass within them is transient, either on its way to 
being accreted or on its way out of the Bondi-Hoyle radius, 
not accumulating over a long term. Our estimates of the ac- 
cretion rate using just the mass that falls onto the particles are 
more accurate. 
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4. THE ACCRETION RATE IN TURBULENT FLOWS 

Thus far we have shown that our model very accurately pre- 
dicts the distribution of accretion rates in terms of the density 
and velocity field of the gas. We can use this result to derive 
an estimate for the accretion rate in a general turbulent flow 
in terms of a few simple parameters of the turbulent region 
that can usually be determined from observations. We con- 
sider a turbulent region of characteristic size I, mean density 
~p, sound speed c s , and three-dimensional Mach number Mo, 
accreting onto an object of mass M. We derive a theoretical 
estimate for the mean and median accretion rate based on the 
probability distributio ns of density, velocity, and vorticity in a 
turb ulent medium in § 14.11 and we compare it to simulations 
in § 14.21 We then discuss the issue of d ifferences between 
mean and median accretion rates in § 14.31 In this Section we 
limit ourselves to discussing accretors randomly sampling tur- 
bulent regions. We discuss how to a pply this results in real 
astrophysical situations in § !5.2l and § 15.31 

4.1. Theoretical Calculation 

Although density, velocity, and vorticity in a turbulent 
medium are certainly correlated to some extent, that correla- 
tion is relatively weak. We can derive reasonable estimates for 
the functional dependence of the accretion rate on the prop- 
erties of a turbulent region simply by neglecting correlations 
and assuming that density, velocity, and vorticity are selected 
independently from the appropriate PDF at every point. 

Numerous authors have stud ied the PDF of densities in a 
turbulent medium dVazauez-SemadenH 119941 iPadoan et alJ 
1997t IScaloetalJ Il998| iPassot & Vazauez-Semadenil 
19981 iNordlund & Padoanl 1 19991 lOstriker et alJ Il999t 
Padoa n & Nordlundll2002l) . IPadoan & Nordlundl (120021) find 
that the PDF is well-fit by the functional form 



dp p (x) 
d\nx 



27TO-2 



= exp 



(lnjc-lnx) 



(16) 



where x = p/~p is the density normalized to the mean density. 
The mean of the log of density (and also the log of the median 
density) is 



hue = 



and the dispersion of the PDF is approximately 



In 1 + 



Mr 



1/2 



(17) 



(18) 



The PDF of velocities is considerably less well-studied in the 
li terat ure. Measuring from the simulations we describe in 
§ 14.21 we find that the PDF is well-fit by a power law with 
an exponential cutoff, 



dp v (u) 
d\nu 



oc w 3 exp(-1.5M' 



(19) 



where u = v/{Moc s ) is the velocity normalized to the velocity 
dispersion of the region. Figure [6] shows our measurements 
from the simulations and our fit. The median velocity for this 
PDF is 1.2Aloc s , but the factor of 1.2 is relatively uncertain 
because it depends strongly on the form of the exponential 
cutoff, which is poorly constrained in our fit. We therefore 
simply take the median velocity to be median(v) = <fi v Moc s , 
where <j> v is a constant near unity. 
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FIG. 6. — Measured values and fit for dp v (u)/d\nu. We show measure- 
ments from our 6 N = 256 runs at A4q = 3 (diamonds), Mq = 5 (asterisks), 
and Mq = 10 (crosses), from our N = 512 run at Mq = 5 (triangles), and our 
best fit (line). 

The vorticity PDF is considerably more complex, as shown 
in Figure In making the Figure we have introduced the 
normalized vorticity 



(20) 



Mac, ' 



which is related to the vorticity parameter for accretion onto 
an object by cj* = uj{Mor^/ '£). The vorticity PDF shows sub- 
stantial variation from run to run, the shape does not appear to 
be independent of Mach number, and the shapes are not easily 
fit by a simple functional form. However, we can still obtain 
useful information by noting that the PDF peaks at roughly 
the same location independent of Mach number, and the me- 
dian normalized vorticity is roughly constant. We estimate it 
as median(cj) = 1O0 W , where </> w is a constant of order unity. 
Note that the fact that this peak is roughly independent of 
Mach number means that the median vorticity in a superson- 
ically turbulent region is completely specified by its velocity 
dispersion and physical size. 

From these PDFs, we can estimate the median accretion 
rate. Since we only expect the fit to be approximate, we sim- 
plify the results by setting A = 1 in equation Q and by droping 
terms of order Mq 2 or smaller when added to terms of order 
unity. With these approximations, and neglecting correlations 
between density, velocity, and vorticity so that we can treat 
the median operator as linear, we predict that the median ac- 
cretion rate, normalized to our simple Bondi-Hoyle estimate, 
will be 



6 med = median(M)/M 
w median 



-1/2' 



/Mo 



1 + Ml/4 



10 / lO^rg 

<f> 6 v V IMP 



-i -0.5 



(21) 
(22) 

(23) 



The functional form agrees with what one would intuitively 
expect. When r^/l is sufficiently small, vorticity has no ef- 
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TABLE 1 

Selected computed accretion rates. 



000.00 



FIG. 7. — Measured values and estimated median for dpu&)/dhiu. We 
show measurements from our 6 N = 256 runs at Mo = 3 (diamonds), Mo = 5 
(asterisks), and Mo = 10 (crosses), and from our N = 512 run at Mo = 5 
(triangles). 

feet on the accretion rate and all vortieity-related terms dis- 
appear. Instead, one is left with a median accretion rate that 
is close to the simple Bondi-Hoyle estimate Mo- It is reduced 
relative to this, however, by the factor (1 + A^o/4)" 1 / 2 , which 
is the ratio of the median to the mean density. When r%/£ is 
large enough that the term in square brackets is much greater 
than unity, the median accretion rate is controlled entirely by 
vorticity. Since the vorticity accretion rate depends approx- 
imately on w~ 9 , and oc <p LJ rs/£, we expect the accretion 
rate scale as (</> u ?"b/^)~ ' 9 , which is exactly what our estimate 
gives. The crossover between the two regimes occurs roughly 
when Mbh = M u , which occurs when the second term in the 
square brackets is about unity. 

We can perform a similar procedure for the mean, although 
here we are considerably hampered by our uncertainties about 
the exact functional form of the vorticity and its dependence 
on Mach number. This uncertainty is amplified by the fact 
that the low vorticity tail of the PDF contributes significantly 
to the mean accretion rate. For this reason, we concentrate on 
estimating the mean a ccret ion rate in the case where vorticity 
is irrelevant, and in § 14.21 we perform a purely empirical fit 
in the case where vorticity is important. For the case where 
vorticity is irrelevant, we compute the mean accretion rate by 
integrating over the density and velocity PDFs, again assum- 
ing that they are independent and setting A = 1 in Q. This 
gives 



:<M)/M 



1 



dp v (u) 



{M, 



2 + u 2\ 3 / 2 du 



du. 



(24) 
(25) 



Although we can evaluate this expression numerically for our 
best-fit PDF, it is more convenient to approximate the inte- 
gral in the case A4o ^> 1 by treating the exponential cutoff 
in dp v (u)/du as a sharp truncation that gives unity at u < <j) u 
and zero at u > </> u , with <fi u ss 1. The factor (j> u represents our 
uncertainty in the exact shape of the exponential cutoff on the 
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velocity PDF. With these approximations, we find 



TtIWA)-!]: 



(26) 



For (f> u = 1, this expression agrees with J25I to better than 10% 
for all Mach numbers between 2 and 200. 

4.2. Comparison to Simulations 

We now compare our theoretically predicted mean and me- 
dian accretion rates to simulations. We run a series of sim- 
ulations of periodic boxes, using the same methodology as 
described in § |3j and compute the mean and median values 
of the accretion rate using our simple model. To get a sense 
of the range of variation in accretion rates, we run six simula- 
tions at a resolution ofN = 256, each using a different random 
realization of the initial velocity field. We examine each run 
at the time when the Mach number has decayed to Aio = 10, 
A^o = 5, and finally .Mo = 3. To compute the accretion rate 
from our simple model, we must also specify r% /£, the size of 
the particle's Bondi radius relative to the size of the turbulent 
region, since the vorticity parameter lj* is proportional to r%. 
We compute mean and cf> me ^ for each realization, at each Mach 
number, and a range of values of / £■ Since our theory will 
probably fail for objects with Bondi-Hoyle ra dii co mparable 
to the size of the entire turbulent region (see § 15. 2> . we limit 

ourselves to r BH /^~0.1, or r B /£< 10 at Mo = 10. 

We report the mean and standard deviation of the six runs, 
for selected values of r^/£, in Tabled an d plot the results in 
Figure[8] The Figure also shows fits based on our theoretical 
predictions for mec | and mea n- F° r me d, we find a good fit 
with equation d23l > for cf> v = 0.93 and (p^ = 1.25. For (/> me an, 
we only have a theoretical prediction for the low rg /£ case, 
so we fit to a function of the form <fi me!ln = ((Mbh)/-^o)[1 + 
f(r%/£, Aio)] q , where f(r% /£,A4o) is a power-law function 
of r B /£ and A4q- We find that the function 



= -|[ln(2<^Mo)- 



1] 1 + 100 



r B 

£M< 



-0.68 



(27) 



with 4> u = 0.95 fits the data reasonably well. Overall, we find 
that our predicted functional forms provide a good fit to the 
simulation data and capture the essential physics of the accre- 
tion process. 

Since the distribution of accretion rates is approximately 



lognormal, we can turn our fits for . 



and me d into fits for 
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FIG. 8. — Measured and fit values of tpmws (upper panel) and me d (lower 
panel) vs. r^/L The plots show runs Mach 3 (diamonds), Mach 5 (asterisks), 
and Mach 10 (triangles), and our fits for these Mach numbers (solid lines). 
The error bars in the upper panel show the standard deviations 111 0mean ■ 

We 

do not show the error bars in me d in the lower panel to avoid cluttering the 
plot. 

the entire PDF of accretion rates. A lognormal distribution 
is completely characterized by two parameters, the center of 
the distribution and its width, and we can solve for these in 
terms of <p mean and </> me d- Doing so, we find that the probability 
distribution of accretion rates in a turbulent medium is 



dp 
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where 
and 



ln(M/M o ) = ln0 med 



<Jm = \ 2\n- 



^med 



(28) 

(29) 
(30) 



4.3. Median versus Mean Accretion Rates 

Since the volumetric median and mean accretion rates we 
have found can be quite different, we would like to determine 
whether one should use a median or a mean accretion rate in 



attempting to follow the evolution of an individual object. The 
mean and median are different because rare, high accretion 
regions contribute significantly to the overall accretion rate, 
even though it is very unlikely for a randomly placed accretor 
to be in one. When an object begins accreting in a turbulent 
medium, its accretion rate is most likely to be near the median. 
Over time, however, it will sample more and more of the tur- 
bulent medium, and its time-averaged accretion rate should 
approach the volume average for the region. We wish to de- 
termine how long this will take, since it is possible that it may 
be longer than the typical time that the object will spend in the 
turbulent region. This will tell us if our volumetric mean and 
median are truly good approximations for the mean and me- 
dian accretion rate in time, measured for a single object rather 
than an ensemble. 

To answer to the question of how long an object must ac- 
crete before reaching the volume mean accretion rate, we con- 
sider the related question of how far an object must move 
through a turbulent volume of gas before reaching the mean 
accretion rate. We answer this using the density and velocity 
field of our N = 512 run at the point when A4q = 5.0, just be- 
fore we insert the sink particles. The rate at which an accret- 
ing particle moves through the volume, accreting from new 
regions, should be roughly equal to the velocity dispersion of 
the gas. Thus, we place a particle at a random point within 
our turbulent box, moving in a random direction with velocity 
v such that |v| = 5.0, while the density and velocity of the gas 
remain fixed. We then compute the time-averaged accretion 
rate along the particle's trajectory, 



M mb (x + \t)dt, 



(31) 



as a function of time T. We repeat this procedure 1000 times. 
Figure |9] shows the fraction of particles for which the time- 
averaged accretion rate is closer to the volume-averaged mean 
than to the median as a function of time. (We use closer in a 
logarithmic sense, i.e. m > v^mean^Wd-) The Figure also 
shows the fraction of particles whose time-averaged accretion 
rates are within a factor of 1 . 1 and within a factor of 2.0 of the 
mean accretion rate in the volume. A majority of the particles 
have time-averaged accretion rates closer to the mean than the 
median after only 0.25 crossing times, where f cr = £/(A4qc s ) = 
0.4, but 95% do not become closer to the mean than to the 
median until more than 6 crossing times. Similarly, more than 
half of the particles are within a factor of 2.0 of the mean 
accretion rate after only 0.28 crossing times, but even after 10 
crossing times fewer than half have time-averaged accretion 
rates within 10% of the mean. Our conclusion is that one 
should use the mean rather than the median accretion rate to 
follow the evolution of any object that accretes for more than 
a quarter of a crossing time, but that one should not expect the 
time-averaged accretion rate to be closer than tens of percent 
to that value unless the particle accretes for a very long time. 

5. DISCUSSION 

5.1. Resolution and Convergence 

It is critical to establish that our simulations have suffi- 
cient resolution for us to believe our results. We first address 
whether we have enough resolution to compute the correct 
accretion rates onto our particles. Our simulation resolves 
the Bondi-Hoyle radius of the particl es by 64 cells, an d their 
Bondi radius by almost 1700 cells. Krum holz et alJ (2004) 
found using the same code that we use here that in simula- 
tions of Bondi-Hoyle accretion a resolution of Ax = tbh/50 is 
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FIG. 9. — Fraction of particles vs. time. The curves show the fraction closer 
to the mean than the median (solid line), within a factor of 2.0 of the mean 
accretion rate (dashed line), and within a factor of 1 . 1 of the mean accretion 
rate (dot-dashed line). 

more than adequate to reproduce the value found by previous 
theoretical and numerical work. Our resolution meets this re- 
quirement, although it is possible that for regions of the flow 
where the velocity is significantly higher than the mean we 
might not have sufficient resolution. At most this will affect 
the low accretion rate tail of the distribution, and should there- 
fore not substantially bias our comparison between the accre- 
tion rate onto particles and our theoretical model. Similarly, 
Krum holz et alJ J2005I) found that Ax = re / 160 was sufficient 
to correctly model accretion with vorticity parameters up to 
= 30. The typical vorticity parameter we find in our simula- 
tion is somewhat greater than this, but by less than an order of 
magnitude. In contrast, our resolution is Ax w rR /1700, more 
than a n order of magnitude greater than that Krumholz et al. 
(2005) found sufficient. 

A related question is whether we have enough space be- 
tween our particles for them to accrete independently of one 
another. A priori we do not expect the particles to interfere 
with one another, since they are separated by 32rBH, much 
larger than any particle's range of influence. There could po- 
tentially be interference between particles if one particle were 
upstream of another and significantly altered the flow of gas to 
its downstream neighbor. However, given the random align- 
ment ot the turbulent flow, having one particle immediately 
upstream of another is extremeley unlikely. Even if one were 
upstream of another, we run the simulation for a time <C 32fBH 
once we insert the particles, so there is no time for depletion to 
propogate from one particle to another. Finally, we can find 
evidence against the possibility of interference between the 
particles by examining the accretion rates we actually mea- 
sure. For Bondi-Hoyle accretion, the accretion rate depends 

3/2 

on the accretion radius asMoc r a c C , and a particle accreting 
at rate Mo has an accretion radius r acc tbh- We would ex- 
pect the particles to begin interfering with one another when 
the accretion radius is comparable to the half the inter-particle 
separation, or 16tbh- This corresponds to an accretion rate of 
64Mq. If interference between particles were a problem, we 




M / M 

FIG. 10. — PDF of accretion rates for the N = 256 runs (dashed lines) and 
the N = 5 12 run (solid line). 

would expect to see the distribution of accretion rates trun- 
cated around 64M(>. However, Figure |4] shows that there are 
no particles accreting at anywhere near this rate, and that even 
our theoretical distribution predicts essentially no particles ac- 
creting at such high rates. We therefore conclude that none of 
our particles have accretion radii large enough for them to in- 
terfere with one another. 

The final question of resolution is whether, before intro- 
ducing our particles, we have enough resolution to model the 
turbulent flow field correctly. This question applies both to 
the N = 512 run into which we insert our particles, and to the 
N = 256 runs we use to estimate the mean and median accre- 
tion rates. To check convergence, we compare the distribu- 
tion of accretion rates we find from the N = 256 simulations 
at Mo = 5.0 to the N = 512 run. We do this for r B /£ = 0.20, 
the value for our particles. Figure Hoi shows the PDFs. The 
N = 512 run mostly falls within the range of variation of the 
N = 256 runs, although it is somewhat steeper, with less of 
the volume at high accretion rates. The average of the mean 
accretion rates in the N = 256 runs is 1 .26 Mo, and the average 
of the medians is 0.36Mo- In comparison, for the N = 512 run, 
the mean is 0.90Mo, and the median is 0.34Mo- The agree- 
ment on the medians is very good, and well within the varia- 
tion among the N = 256 runs, but the agreement on the means 
is only fair. The difference in quality of convergence proba- 
bly occurs because the medians are sensitive to the peak of the 
distribution, but the means depend more on the low vorticity, 
high accretion rate tail of the distribution that is more sensi- 
tive to resolution. Thus, our means probably have a ~ 25% 
systematic error. 

As a final note, the behavior of vorticity with resolution is 
interesting, and helps explain why we are better converged in 
the median than in the mean accretion rate and why the mean 
is decreasing with resolution. We compute the mean and me- 
dian vorticity in the velocity field of our N = 512 run when 
.Mo = 5.0, just before we insert the sink particles. We then 
sample the simulation down to N = 256, N = 128, N = 64, 
and N = 32 and compute the mean and median vorticity for 
those resolutions. Both the mean and median vorticity in- 
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FIG. 1 1 . — Mean (upper tine) and median (lower line) normalized vorticity 
ui versus resolution. 

crease monotonically with N, but the median vorticity appears 
to be converging while the mean only shows slight curvature 
and may actually be diverging. We show this in Figure ITU 
The mean is controlled by a high vorticity tail, and it appears 
likely that there is no converged value for it until one reaches 
the dissipation scale, either due to grid viscosity (in a simu- 
lation) or due to physical viscosity (in a real fluid). In con- 
trast, the median appears to be well-defined and converged at 
N = 512; the difference in median between the N = 256 and 
N = 512 resolutions is less than 2%, and the median vorticity 
is rising with N only as N - 026 . A similar trend is visible in 
Figure where the N = 512 run peaks at roughly the same 
location as the N = 256 runs, but has a longer high-vorticity 
tail. 

The divergence of the mean does not affect the convergence 
of our results, however. The high vorticity tail that causes the 
mean vorticity to diverge corresponds to low accretion rates. 
Thus, it will not cause the mean accretion rate to diverge. 
Nor will it cause our accretion rate to go to zero as resolu- 
tion increases. The high vorticity tail will simply contribute 
less and less to the mean accretion rate as the resolution in- 
creases, but the contribution from vorticities near the median 
will not decrease because the median does not change with 
resolution. Thus, we can regard our ~ 25% change between 
the N = 256 and N=512 runs as a rough estimate of the likely 
error from imperfect convergence. However, we conjecture 
that vorticity on scales much smaller than tbh will not sig- 
nificantly affect the accretion rate. In our N = 512 runs, we 
have tbh/Ax = 4 on the base AMR level. (Once we insert the 
particles we greatly increase the resolution around them, but 
our estimates of </> mea n are based on the simulation cube before 
that insertion). Since we already have marginal resolution at 
N = 512, we believe it is unlikely that increasing resolution 
past N = 512 will stronlgy affect me an- Also note that this 
convergence issue only affects accretion when /£ is large 
enough to affect the accretion rate. The velocity PDFs, as 
shown in Figure[6] are already very well converged at N = 256, 
so only when vorticity becomes important are our results un- 
certain. 



5.2. Limitations 

There are several limitations to our results. First, we have 
neglected magnetic fields. These are likely to suppress the ac- 
cretion rate relative to what we have found, because they will 
restrict the flow of gas across magnetic field lines to an ac- 
creting object. In molecular cloud s, where the Alfven Mach 
number is of o rder un ity (lMcKedll989t [ McKee et alJll993l 
ICrutcherll99^) or less iPadoan et alJ2004l) . the effect is likely 
to be of order unity. However, magnetohydrodynamic simu- 
lations similar to the ones we have performed will be required 
to make a more definitive statement. A related limitation is 
that we have not considered a self-gravitating medium. If the 
gas is strongly self-gravitating, its density and velocity distri- 
butions will change, and our results for mean and ^> me( j may 
need to be modified. However, our central result that vorticity 
is critical in setting the accretion rate will not be affected by 
the addition of self-gravity. 

A second limitation is that our approach is likely to break 
down for subsonic turbulence. The distribution of densities 
and velocities is different in subsonically turbulent media than 
in supersonically turbulent ones, so our results will not ex- 
trapolate to that regime. For very subsonic turbulence where 
there is some vorticity present, it is likely that the accretion 
flow approache s the vorticity-d ominated low cj* regime iden- 
tified by Krumholz et al. (2005). The transition should occur 
for values of A^o near unity. However, we have not mapped 
out the details of the change from one regime to the other, so 
we must be careful in attempting to apply our theory to the 
warm atomic or ionized phases of the ISM, since the motions 
in these phases are generally transsonic or subsonic. 

Third, our results are only applicable to objects with Bondi- 
Hoyle radii significantly smaller than the size of the turbulent 
gas cloud in which they are accreting. Our estimate of the 
accretion rate is based on the assumption that, on large scales, 
the turbulent gas is uniform and can be roughly modeled by a 

periodic box. However, if 7bh ^ £, then the turbulent gas will 
become centrally condensed due the gravity of the accretor, 
and our assumption of uniformity will fail. 

The final major limitation is that we have only considered 
the case in which the accreting object is not moving through 
the turbulent medium at a velocity comparable to or larger 
than the velocity dispersion of the gas. In such a situation an 
accretor should undergo velocity-dominated accretion regard- 
less of its position in the turbulent medium. The accretion rate 
will be the standard Bondi-Hoyle rate, using the local den- 
sity as poo. Consequently, the distribution of accretion rates 
should be proportional to the distribution of dens ities, which 
is we ll measured from simulations (e.g. IPadoan & Nordlundl 
1999). This limitation prevents us from applying our theory to 
neutron stars or black holes passing through molecular clouds, 
since these objects generally have velocity much larger than 
the turbulent velocity dispersions of GMCs. 

One might think that our assumption of isothermality also 
presents a considerable limitation. However, this is not likely 
to be the case. In ordinary Bondi-Hoyle accretion, changing 
the equation of s tate only changes the accretion rate b y factors 
of order unity (lRufferll994lRuffert & Arneti 19941) . Thus, a 
different equation of state is not likely to change the nature of 
the accretion flow onto a particle. The only way that a change 
in the equation of state might affect our accretion rates is by 
changing the distribution of densities, velocities, and vortici- 
ties in the turbulent gas. Since we are assuming the turbulence 
is highly supersonic, a different equation of state is unlikely 
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to affect the velocity field very much; thermal energy will al- 
ways be considerably less than kinetic energy over most of 
the flow. The remaining potential effect of non-isothe r mahty 
is on the density field. Passot & Vazquez-Semadeni ( 1998) 
show that isothermal turbulence leads to a PDF of densities 
that is lognormal. Polytropic indices 7 < 1 lead the PDF to 
develop a power-law tail at high densities, and indices 7 > 1 
lead to a power-law tail at low densities. The latter is unlikely 
to affect the accretion rate, since the low density gas that is 
affected does not contribute significantly to the mean or me- 
dian accretion rate. If 7 < 1, there will be more gas at high 
densities than in the isothermal case, and the accretion rate is 
likely to be larger than we have estimated by an amount that 
depends on how small 7 is. Our results are therefore likely to 
change only by order unity for 7 > 1 , and will only change by 
much more than that if 7 <C 1 . 

5.3. Implications 

Our results have implications for several astrophysical 
problems. Here we sketch out t hree. The firs t is the problem 
of pre-main sequence accretion. Padoan et al. ( 2005 ) consider 
accretion onto T Tauri stars from a turbulent medium using the 
ansatz M tm b = Mbh- (Note that this is slow accretion that takes 
place after the pre-stellar core has collapsed and been ac- 
creted, not the competitive accretion model of star formation. 
We discuss competitive accretion below.) They find an ac- 
cretion rate consistent with observed accretion rates onto pre- 
main sequence stars and brown dwarfs in nearby low-mass 
star-forming regions. We have shown that using M tur b = Mbh 
can lead to significant overestimates of the accretion rate, and 
that the mean and median accretion rate are quite different. 
Accretion rates are inferred from observations of the accre- 
tion luminosity, which is sensitive to the accretion rate on time 
scales far smaller than the crossing time of the system. Thus, 
one should use the median rather than the mean accretion rate 
when comparing to individual accreting objects. However, 
we find that Padoan et al.'s conclusions are not affected too 
strongly by our results, provided their application is limited 
to the low-mass star-forming regions to which Padoan et al. 
compare their models. Such regions have widely distributed 
pre-main sequence stars and star-forming clumps with large 
radii. If stars are forming in a turbulent clump of size I ~ 5 
pc, as Padoan et al. assume, then for a 1 Mq star r% /£ ss 0.02 
(for 10 K gas), and vorticity is relatively unimportant. 

On the other hand, massive star-forming regions are more 
compact and turbulent, and star-forming clumps there are 
more typi cally on scales £ s a 0.5 pc, with Mach numbers 
Mq ~ 20 (Plu me et allll997l) . Levels of turbulence are not 
significantly lower in starless cores, so it is likely that these 
Mach numbers ar e typical even at early stages of the star 
formation process ( Yonekuraet al. 2005). For such a region 
r^/lK, 0.2 for a 1 M star, and vorticity can suppress accre- 
tion. This conclusion leads to the possibility of an interesting 
observational test of our conclusions. Assuming that Padoan 
et al. are correct and that accretion onto protostars comes 
primarily from accretion of gas from the environment (rather 
than from some internal process involving an accretion disk), 
one could bin accreting stars by mass, and look for a change 
in the dispersion of accretion rates versus mass. This break 
should occur because the mean and median accretion rate be- 
gin to be affected by vorticity at different values of r% / 1. Us- 
ing J27i . in a clump with Ado = 20, 4> me:m is independent of 
mass for stars of mass M < 1 Mq, and </> me an ^ M -0 - 7 for more 
massive stars. On the other hand, </> me d is independent of mass 



up to a mass of 130 M Q , vastly larger than any observed pro- 
tostar. Thus, above 1 M Q the ratio </> me an/Vmed oc M" 7 , and 
the dispersion of accretion rates, which depends on this ratio, 
should decrease with mass. Given a reasonably large sam- 
ple of accreting protostars in a dense clump with a range of 
masses above and below 1 M , one could see this decrease 
in dispersion with mass. One could also in principle look for 
the change in cf> meim with mass directly, but this is more diffi- 
cult because one would need a very large sample of accreting 
stars to determine mea n accurately. Measuring the dispersion 
requires fewer sources. Even so, measuring accretion lumi- 
nosities of T Tauri stars in massive star forming regions is a 
significant observational challenge due to distance, confusion, 
and extinction. Unfortunately one cannot test our model using 
more massive stars, because the rate of Bondi-Hoyle accretion 
onto stars larger than ~ 10 M^ is likely to be subst antially re- 
duced by radiaton pressure (Edgar & Clarke 2004). 

A second potential implication is in the area of competi- 
tive accretion. In this model of star formation, the initial 
mass function is determined by seed protostars accreting un- 
bound gas within the turbule nt molecular clump where they 
form ( e.g. Bonn ell et al ll997tlKlessen & Burkertl2000l.l2001t 
iBonnell et al. 2001a b, 2004). Our theory does not apply dur- 
ing the first stage of competitive accretion, when a star is still 
accreting from its initial bound core, for two reasons. First, 
the core is bound to the star and to itself, while we have 
considered accretion of unbound gas. Second, the star is co- 
moving with the gas in its core, so it is not randomly sampling 
the turbulent flow. Once the star accretes its parent core, these 
limits no longer apply and one may use our theory. The star is 
no longer accreting bound gas, and, since the time required to 
accrete the bound core (which is comparable to the dynami- 
cal time of the parent clump - McKee & Tan 2003) will allow 
the star to become dynamically decoupled from the turbulent 
flow, the star will be randomly sampling the gas in the molec- 
ular clump to which it is gravitationally confined. 

Our results suggest that in this stage the accretion rate in 
such clumps may be considerably different than a naive ap- 
plication of the Bondi-Hoyle formula suggests. Simulations 
of competitive accretion performed to date likely have too lit- 
tle resolution to accurately model the effect we describe here. 
In a typical massiv e star-forming clump, the velocity disper- 
sion is ~ 4 km s" 1 dPlurne et al.lll997tlYonekura et al1 l2005). 
giving a Bondi-Hoyle radius for a 0.1 M seed protostar of 
only 5.5 AU, far smaller than the resolution of the simula- 
tions performed thus far. Published calculations of competi- 
tive accretion to date have largely avoided the problem by us- 
ing the lower velocity dispersions characteristic of low-mass 
star forming regions. However, it represents a formidable ob- 
stacle to simulating the formation of rich clusters from mas- 
sive clumps. It is unclear what effect failing to resolve this ra- 
dius would have on the computed accretion rate. Competitive 
accretion simulations generally use smoothed-particle hydro- 
dynamics with sink particles, and these methods compare the 
kinetic and potential energies of a gas particle in an attempt to 
determine whether it is bound to a sink particle before accret- 
ing it (Bate et al. 1995). As a result, the most likely error is an 
understimate of the accretion rate caused by failure to resolve 
the Bondi-Hoyle accretion shock cone and the resulting loss 
of gas kinetic energy. 

A third potential implication of our result is for the accre- 
tion rates of black holes in the centers of galaxies. Some sim- 
ulations and semi-analytic models of the growth of galactic 
center BHs assume that the black holes accrete at the Bondi 
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rate (e.g. ISpringel et aLll2005l) . If the gas around the black 
hole is turbulent, however, using the Bondi rate may lead to 
significant errors. It is uncertain whether the accretion process 
onto galactic center BHs can be described by Bondi accre- 
tion at all, but, if it can, then simulations and analytic models 
should use the turbulent accretion rate we have found rather 
than the unmodified Bondi rate. 

6. CONCLUSIONS 

The primary result of our investigation is that the accretion 
rate in a turbulent medium is significantly enhanced over the 
zeroth-order estimate one obtains by simply using the turbu- 
lent velocity dispersion in the Bondi-Hoyle formula. How- 
ever, this enhancement may be partially or completely offset 
by the effects of vorticity within the turbulent medium, and 
may even be suppressed below the Bondi-Hoyle rate for Mach 
numbers near or below unity. For this reason, simply using 
the Bondi-Hoyle formula either by treating the global veloc- 
ity dispersion as a sound speed or by appl ying the fo rmula 
to every point within the medium (Padoanet al. 2005j) does 
not produce the correct accretion rate. On the other hand, 
approximating the accretion process as being dominated by 
either velocity or vorticity, and using the lesser of the two im- 
plied accretion rates, appears very consistent with the results 
of simulations. 

We use this simple model to compute the probability dis- 
tribution of accretion rates onto point particles in a turbulent 
medium. The PDF of accretion rates is lognormal, and we 
give simple formulae that specify the shape of the lognormal 
distribution, and the mean and median accretion rates, over a 
relatively wide range of accretor masses and turbulent Mach 
numbers. The mean and median accretion rates can be quite 
different, due to the long tail of the lognormal accretion rate 
distribution, and it takes an accreting particle of order a cross- 
ing time of the turbulent region before its time-averaged ac- 
cretion rate approaches the mean for the region. 

Our findings suggest that a number of previously published 



numerical and analytical results may need to be reconsidered, 
and that future work may have to meet quite severe resolu- 
tion constraints. Estimates of accretion rates in star-forming 
regions using the simple Bondi-Hoyle formula are probably 
incorrect, and simulations of competitive accretion that do not 
resolve the Bondi-Hoyle radius, and possibly even smaller 
scales, are likely to produce incorrect accretion rates. Sim- 
ulations of star formation that have been completed to date 
have largely avoided this problem by not simulating the very 
large velocity dispersions found in massive star-forming re- 
gions, and our work shows that expanding the simulations to 
these regions will require much higher resolution. We leave 
a more detailed reconsideration of these problems, using our 
newly developed theory of Bondi-Hoyle accretion in a turbu- 
lent medium, to future work. 
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